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Abstract 



Parametric resonance has been discussed as a mechanism for copious par- 
ticle production following inflation. Here we present a simple and intuitive 
calculational method for estimating the efficiency of parametric amplification 
as a function of parameters. This is important for determining whether reso- 
nant amplification plays an important role in the reheating process. We find 
that significant amplification occurs only for a limited range of couplings and 
interactions. 
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Reheating after the end of inflation is the process during which almost all elementary 
particles in the universe are presumed to have been created. In versions of inflationary 
cosmology in which inflation ends by slowly rolling towards the minimum of the inflaton 
potential,0'0 reheating occurs as the inflaton field descends to the minimum and begins to 
oscillate rapidly about it. The oscillating inflaton field decays to the various matter fields 
to which it is coupled.&^i In the typical case, the decay of the inflaton can be represented by 
adding to the equation-of-motion for a phenomenological term Tcf) where V is the effective 
decay rate of particles; a corresponding term is added to the equation describing radiation 
to account for the transfer of energy from to radiation.! This treatment assumes the 
particles composing the oscillating 0-field decay incoherently. 

Recently, there has been interest in a two-stage reheating scenario which begins with 
the coherent transfer of energy from the oscillating inflaton to a second field (or fields) 
X through parametric resonance. Later, the X-particles decay to a rapidly thermalizing 
distribution of ordinary matter, completing the reheating process. During the first stage, 
which has been named preheating, the density of X particles undergoes a period of exponen- 
tial increase. The resonant transfer of energy is important to inflationary cosmology if it is 
"efficient" : if it results in amplified production of X particles by many orders of magnitude 
compared to estimates based on the decay rate (T) of an individual inflaton particle such 
that a large fraction of the inflaton energy density is converted to X particles. In addition 
to producing a higher reheating temperature, preheating enables the coherent transfer of 
inflaton energy to particles X which are more massive than itself, whereas only particles 
less massive than are produced in incoherent decay and the density of very massive X- 
particles is exponentially suppressed during subsequent thermalization. This makes possible 
models of baryogenesis which rely on there being a high density of particles more massive 
than following inflation which subsequently undergo baryon-asymmetric decays.0 

The goal of this paper is to present a straightforward and intuitive method for de- 
termining the efficiency of parametric resonance for a wide range of model couplings and 
interactions. We present results for a simple model in which the X-particle is massive and 
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has quartic self-interactions. The method predicts and makes intuitively clear why resonant 
particle production is efficient only for a narrow range of masses and couplings. To check 
the prediction, we have developed precise numerical codes similar to those of Tkachev and 
Klebnikov0El and, as shown in the paper, we have found that the estimates of particle pro- 
duction compare very well. Where there is overlap, the results are also consistent with the 
analysis by Kofman et al. .E3IH1 With this verification, this approach can now be applied to 
general inflaton models to determine if resonant particle production is a significant feature 
of the reheating process. 

In our method, resonance is described by trajectories through stability/instability regions 
of the Mathieu equation. This notion has been alluded to in previous discussions JllllQ' III 
here we show that it can be developed into a reliable, predictive tool. Quantitative estimates 
are simply obtained from calculable properties of the Mathieu equation and the trajectories. 
The estimates are reliable except for very large couplings such that other nonlinear effects - 
stochastic resonance, backreaction, rescattering and chaotic instability — become important. 
Although this large coupling regime has been the focus of many interesting studies,B000'lll 
it lies far into the "efficient preheating" regime; as we will demonstrate, we do not need to 
consider these complications for our purpose of finding the boundary between inefficient and 
efficient preheating. 

The general model which we consider is described by the Lagrangian density: 

£ = + \iW 2 - \^ 2 X 2 - \rny - Ir&C ~ V{X). (1) 

Here V(X) represents self-couplings of the X-field of cubic order or higher. If we assume 
that is spatially uniform, we can derive a set of coupled equations for <fi and the Fourier 
modes , X^. The equations for X^ decouple from one another in the Hartree (mean field) 
approximation. In terms of "comoving" fields, $ = a 3 ^ 2 cf> and Xk = a3 ^ 2 Xk, the linearized 
equations are: 



ml + + 6nGp 



$ = (2) 
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where 
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Xk = 0, (3) 



u 2 k = ^ + ml + f(x)+&KG P (4) 
ct 

1 ,'d la 2 . 

P= -^G { a + 2a 2) (5) 
a 3 / 2 V'(y/a 3 / 2 ) 

fix) = ( W J ) (6) 
X 

where p is the pressure, a is the scale factor, G is Newton's constant, and (. . .) represents 
the spatial average. Henceforth, we choose units of energy, mass and time where = 1. 

The pressure is negligible when the universe is dominated by inflaton oscillations. Fur- 
thermore, until one reaches the condition where g 2 x 2 /a 3 ~ m/i, the backreactioni@0lll'lll 
on $ due to X and rescattering are negligible, too. This condition is not satisfied as long 
as the energy density in X-particles, p x = oj 2 {\ 2 /a 3 ) = G(l)(g 2 $ 2 + m 2 )(x 2 /a 3 ), does not 
exceed the inflaton energy density, m\<& 2 . But, this is precisely our regime of interest since 
we wish to determine the range of parameters for which resonant amplification is inefficient 
to barely efficient (p x < p^). Hence, for our purpose of finding the boundary between the 
inefficient and efficient regimes, backreaction and rescattering can be safely neglected. Of 
course, for previous studies in the highly efficient, strongly coupled regime, backreaction and 
rescattering cannot be ignoredBBBSEiJ 

Hence, to good approximation, $ behaves as a simple harmonic oscillator, $ = <£>ocost, 
where time t is measured in units of = 1 beginning from initial time to ~ m^ 1 when the 
inflaton begins to oscillate. Then, the equation-of-motion for x can De written: 

^ + [A + 2gcos(2t)] Xfc = (7) 

where 
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A = u 2 k+ 9 ^ = u 2 k + 2 q (8) 



4o 3 

"2a 3 
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In the limit of a static universe, a — > constant, Eq. fl7|) is precisely a Mathieu equationS 
The (g, A) plane is separated into stable regions where the amplitude of x stays constant 
and unstable regions where the amplitude increases exponentially, x oc e v< , where the critical 
exponent v (the imaginary part of the Mathieu characteristic exponent) depends on q and 
A. See Figure 1. For fixed A, v grows and the instability strips broaden as (q, A) moves 
horizontally towards the right. The dotted line A = 2q roughly divides the strips into their 
narrow and broad regimes. Each mode k corresponds to some fixed (q, A); if this point 
lies in an unstable regime, the mode is exponentially amplified. Increasing k corresponds to 
increasing A leaving q unchanged, so the modes lie along a ray pointing vertically upwards 
beginning with the k = mode. 

Important differences arise in an expanding universe. First and foremost, A and q are 
functions of the scale factor a and, therefore, vary with time. In our treatment, we propose 
to use the Mathieu equation as a guide but take account of the time-evolution in the (q, A) 
plane. For simplicity, we first consider f(x) = 0. As shown in Figure 1, (q = qo/a~ 3 , A = 
2q + (fc/ct) 2 + m 2 ) is no longer a fixed point in an expanding universe, a = a(t). Rather, 
each k corresponds to a trajectory in the (q, A) plane from upper right, (q , A ), to lower 
left, (q, A) — > (0,m^). For k > 0, the starting point is shifted vertically upwards from 
(A = 2q + m 2 x ) by k 2 . Then, the mode is red shifted and the trajectory curves towards the 
A = 2q + ml line as q = g 2 ^l/Aa 3 — > 0. 

Resonant amplification occurs as the trajectories pass through parts of the instability 
regions with large critical exponent v. The greatest resonant amplification occurs for a small 
band of modes near k = (see Figure 1) since the corresponding trajectories lie furthest 
to the right in the (q, A) plane passing through regions with the highest values of v. Once 
a trajectory passes below some critical value, q = q e , though, the value of v approaches 
zero and there is no further amplification. This means that there is only a finite time for 
resonant amplification in the expanding universe, whereas there is no limit in the static case; 
this is a key difference that qualitatively changes the nature and robustness of parametric 
resonance. For example, a consequence is that the coupling g must exceed a certain minimal 
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value for there to be significant amplification, as we shall show. Another prediction is that 
increasing m 2 . > (or /(%) > 0) suppresses resonant amplification since it moves the band 
of trajectories further to the left towards the more stable regimes. 

Our Mathieu equation picture can be made quantitative and used to estimate the energy 
density in X, p x (t), produced by resonant scattering as a function of parameters. The 
zero-point contribution of mode k to the energy density is Pk(t) = where ftf.(t) = 

u;l(t) + (g 2 $ 2 /a 3 ). The coupling of $ to X appears as an additional effective mass for the 
X field. The production of particles by parametric resonance adds energy density 



beginning from time t — ^i 171 ^ 1 ) — when inflation ends and inflaton oscillations 

begin; here, v^if) is the critical exponent for mode k at time t. 

Only modes within a small band near k = contribute significantly to the final, 
integrated p x (t). See Figure 1. For non-zero comoving wavenumber k, the trajectory 
A = 2q + {k/a) 2 + m 2 , > 2q lies to the left of the A = 2q line where there are smaller 
values of the critical exponent u k {t) on average and, hence, there is less amplification. As 
the trajectory proceeds to smaller q, the value of u k averaged over several instability regions 
(= Ukif)) increases as k/a decreases. To gauge the increase in &k(t), it is useful to consider 
the family of curves of the form A = 2q + \/bq, along each of which the time-averaged 
v = v(b) is uniform {e.g., as can be confirmed using Mathematica S.Otm)- The value of 
9(b) increases as b decreases, reaching up to v(b) ~ 0.1 for 6 = 0. An important divider 
is the curve with b = 0{1) for which V(b) pa 0.05. Modes with trajectories lying above 
this curve have very small critical exponents, V k < u(b) <^ 0.1, and so undergo very little 
amplification. Hence, only the small band of modes close to k = with trajectories lying 
below this curve are significantly amplified. The resonant band corresponds to < k 2 < k 2 
where k 2 = \fhq~o — m 2 ,, which decreases slightly as m x increases. (A similar relation has 
been introduced based on different arguments in Ref. 6.) Also, Uk=o decreases slightly as 
m x increases because the k = mode is shifted to the left in the massive case. We can 
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(9) 



approximate / dt'ukit') pa Vk=o(t — t ) and obtain a simple expression for the total amplified 
energy density: 



where Uk=o ~ 0.1 for m x <C 1 and and Uk=o ~ 0.075 for m x > 1. 

The key parameter, V k=Q , can be estimated without recourse to numerical averaging over 
a trajectory once one is familiar with some basic properties of the Mathieu plot. From the 
Mathieu plot, one can see that v varies between (its value in the stable regimes) and 
0.27 (its maximal value in the unstable regimes). For m x = 0, the trajectory moves through 
stability and instability regimes such that the stability regimes are wider than the instability 
regimes (that is, for the range of small g's responsible for amplification). Hence, one can 
anticipate that the average v has a value somewhat less than half its maximum (0.135). 
Indeed, a numerical integration along the trajectory in the (q, A) plane yields the result 



For the massive case (m x > 1), it is evident that v has to be smaller than the massless 
result, 0.12: the trajectory lies further to the left of the A = 2q line where the critical 
exponent is lower. On the other hand, as shown in Figure 1, the trajectory lies to the right 
of the curve A = 2q + y/bq. Along this curve, v averaged over several stability/instability 
regimes is constant; for b = 1, say, the average value of 0.05. Hence, one can anticipate 
that the average value for m x > 1 lies between 0.05 and 0.12. Numerical integration shows 
v ~ 0.08, in agreement with this estimate. 

Resonance continues until q(t) = qo/a 3 = qo(to/t) 2 falls below the critical value q e . 
For the massless case, q e ~ 1/4 below which trajectories enter a Mathieu stability region 
that sustains as q — > 0. For the massive case (m x > = 1), resonance sustains until A = 
2q e + \Jb~q e = 2q e + ml, or q e = \m x . The duration of the resonance is t — 1 = {\JqoJq~ e — l)^o- 
To reach barely efficient amplification (p x pa p^), it is necessary that go Qe- 

Putting these details together, we obtain our prediction for general m x 
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V = 0.12. 



[(ko) 1/2 




(m 2 x + ^^(toAe) 2 ) 1/2 (exp[0.2 V ^7J* ] - 1), 



(11) 
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where go = 9 2 ^o/^, Qe = max{m^/6, 1}, t e = yqo/q e to and to ~ 1- For the case of 
interactions, /(x) > 0, we replace rr? x — > + f(x)- 

Figure 2 compares our prediction to the results of our exact numerical code based on 
Eqs. (|) and (§). The figure shows an example of efficient preheating, p x (t e ) / Ptotai — ^(1)> 
and several choices of parameters which produce inefficient preheating, p x (t e ) <C p^ < ptotai- 
In both regimes, our heuristic picture predicts the general behavior and estimates p x (t e ) / p to tai 
reasonably well. We include an example with V(X) = |AX 4 , a case of self-interactions 
which has been discussed qualitativelyElill and computed numerically.@@ Here we are able 
to predict quantitatively its effect simply by modifying our expressions for q e and Q fc=0 (t) 
in Eq. ([Tl|) to take account of fix)- 

The numerical solutions show fine and intermediate scale structure which can also be 
understood in our Mathieu picture. Consider, for example, the x particle number ( 

ti^) as a 

function of time, as shown in Figure 4(a). Because n x oc l/fi x oc 1/$, there is a sharp spike 
each time $ oscillates through zero, the minimum of its potential. This is an artifact of the 
definition; the physically important issue is how the particle number changes between spikes 
and from before to after each spike. At early times when q is large, one observes that n x (t) 
is flat between spikes and undergoes a discrete jump up or down (mostly up) from one side 
of the spike to the other. In discussing the strong coupled regime, Kofman, et a/0 have 
pointed to this behavior and argued that the up and down jumps become stochastic when 
q passes through many stability/instability bands within one inflaton oscillation and the 
coupling is large. However, here we see the same jump feature in the more weakly coupled 
regime even though q is not changing as rapidly. And, while there are occasional downward 
jumps, the average behavior is not stochastic; they are well fit by an exponential envelope. 
This is because the jumps and plateaus themselves have nothing to do with stochasticity or 
rapid evolution of g; they are already features of the static (constant q) Mathieu equation 
when q is large. The discrete jumps are illustrated in Figure 4(b) for constant q = 119. 
Because q is constant, all the jumps are upward. In numerical experiments where q varies 
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slowly from a instability regime to a stable regime during an oscillation, we find downward 
jumps. Hence, all the small-scale structure in the Figure 4(a) can be understood in terms 
of the static or near static Mathieu equation. 

Figure 4(c) shows the solution for the static Mathieu equation when q is small (q = 0.6), 
in which case the particle production is found to increase continuously between spikes and 
across spikes. This behavior compares well with the late time behavior in Figure 4(a) when 
q is small. 

This behavior can be understood semi-analytically. In either the high or low q regime, 
the solution of the static Mathieu equation is described by a Floquet solution, 



That is, there is an exponential envelope multiplying a modulating function, /(£). The 
difference in the low- versus high-g limit is that the modulating function is continuously 
increasing at low q whereas it has sharp jumps at high-g. The particle number, then, is 



where g(t) is a function of /, /', /" which has a step-like behavior for high q and is continu- 
ously increasing for low q. All this behavior is reflected in the expanding universe solutions 
found numerically. While the details of the modulating function are interesting and it is 
gratifying that they can be understood within the Mathieu picture, what is important for 
our purpose of predicting the average particle production is the fact that the envelope func- 
tion is exponential and that v is predictable by taking the average over trajectories. Figure 
4 illustrates that the exponential envelope as predicted by the Mathieu picture is a good fit 
in the static and expanding universe limits in both the high and low q regime in the weak 
coupling limit. 

In addition to small-scale structure (on the scale of a single oscillation) in at early 
times, the particle number n x in Figure 4 and the energy density p x in Figure 2 also show 
intermediate scale structure — a sequence of rises and broad plateaus stretching over many 



X = e» t f(t). 



(12) 



n x = e 2 ^g(t) 



(13) 
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oscillations. These occur because g decelerates as the universe expands and spends several 
oscillations in each of the last few stability and instability regions. The rises/plateaus corre- 
late with individual stability/instability bands. During this period, most of the amplification 
takes place. 

We have argued that the Mathieu picture explains qualitatively all the features seen in 
particle production in the inefficient to barely efficient regime. We have also seen from our 
Figures that the picture agrees well quantitatively with numerical results, in spite of the 
fact that it ignores stochasticity, backreaction and rescattering effects which are included 
in the code. The code itself is standard. It employs a variable time step size Runge-Kutta 
algorithm. We approximate the Fourier transform of a field by a discrete lattice in k-space. 
Each lattice point corresponds to one of the set of coupled differential equations. The 
rescattering is taken into account by splitting the inflaton field into a zero-mode piece and a 
fluctuation piece: = <po+5(p. We treat 0o as a classical background field, and decompose Scf) 
into k-modes, with a treatment similar to the of the \ field. The good quantitative agreement 
between the Mathieu picture and numerics found support our theoretical argument that 
backreaction and rescattering are not important for determining the boundary between 
inefficient and efficient preheating. They only become non-negligible for couplings far into 
the efficient preheating regime. 

Our analysis predicts that non- negligible resonance occurs for only a limited range of g. 
First, there is a minimum value of go needed to have non- negligible resonance, which implies 
a lower bound on g. The values of and $o are both fixed by the constraints of having 
sufficient e-folds of inflation and sufficiently small density fluctuations. The only freedom left 
to adjust go = 9 2 &o/4 is through the coupling g. In order to have go > g e , it is necessary that 
g 2 $g > max{|m^., 1}. However, there is also an upper bound on g. The resonance picture 
assumes g 2 $l <C M 2 or else the mass of the X-field is (m 2 + g 2 $ 2 ) 1//2 > M p and quantum 
gravity effects become important. In our example, chaotic inflation for V{4>) = ^m^cp 2 , 
$ ~ M p at the end of inflation, so we must have g <C 1. Also, g < 10~ 3 is required in order 
to avoid large radiative corrections (assuming no supersymmetric cancellations) .0 
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A second prediction is that resonant amplification falls sharply as m x increases above 
the inflaton mass. For example, in our model, inflation requires « 10 13 GeV. A proposed 
application of parametric resonance has been to produce A-bosons with mass 10 14 GeV or 
higher whose decay may produce baryon asymmetry Our analysis suggests that parametric 
resonance is not efficient for m x jm^ 3> 1 unless g is much larger than in the massless case. 

In general, our method can be used to determine the range of g and m x resulting in 
efficient preheating, as shown in Figure 3. Inefficient preheating means that parametric 
resonance is inconsequential and reheating occurs predominantly through incoherent decay. 
This curve shows that the efficient preheating regime is small - the minimal value of g 
needed to have efficient preheating rises sharply with m x , ultimately exceeding g = 1 (where 
all approximations become invalid) at m x rs 100. Hence, parametric resonance cannot be 
relied on to produce A particles many times more massive than the inflaton field,@ which is 
problematic for models of high energy baryon asymmetry generation. Adding interactions 
with /(%) > further suppresses resonance. In particular, efficient reheating requires a 
small quartic coupling, A <C g 2 (see Figure 2). 

Although our results have used explicitly a quadratic inflaton potential, we have checked 
that the trajectory picture works well a quartic inflaton potential as well, in which resonant 
amplification can also be described in terms of stability-instability regions. El In either case, 
an issue that remains to be explored is an apparent, classic chaotic instability that occurs 
in our numerical codes for choices of g and m x far into the efficient preheating regime of 
Figure 3 as backreaction and rescattering become significant. 

In summary, we have developed a picture of parametric resonance in an expanding uni- 
verse in which resonant A-particle production is associated with passage through a sequence 
of instability regions of the Mathieu equation. The picture explains qualitatively and quan- 
titatively why resonance in an expanding universe is effective for only a constrained range of 
couplings, masses and interactions, especially in the limit of large m x . In particular, we have 
shown that efficient preheating occurs only if the dimensionless couplings satisfy somewhat 
unusual conditions where the coupling g is large but all other interactions of the inflaton 
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and X-field are small. Unfortunately, this means that parametric resonance in inflationary 
cosmology is less generic or robust than was hoped for. 

We would like to thank L. Kofman, A. Linde, D. Boyanovsky, I. Tkachev and R. Hol- 
man for useful discussions in the course of this work. This research was supported by the 
Department of Energy at Penn, DE-FG02-95ER40893. 
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FIGURES 



A 




FIG. 1. A sketch of the (q, A) plane in which grey (white) strips represent instability (stability) 
regions of the Mathieu equation. The dotted line corresponds to A = 2q. The thick solid line 
corresponds to A = 2q + m^, the trajectory of the k = mode (assuming V(X) = 0). Also 
sketched (see thin solid lines with arrows) are trajectories with k < k, which dominate the resonant 
amplification. Only trajectories lying below the dashed curve, A = 2q + ^/bq with b = 0(1), are 
significantly amplified; in the text, we explain how to use this curve to determine k and q e . 
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FIG. 2. Plot of Px/ Ptotai vs. m^t/27r comparing numerical results (thin, ragged curves) to the 
heuristic Mathieu equation picture for various models (thick, smooth curves), where 27r/m^ is the 
inflaton oscillation period. The approximation reliably predicts the time-averaged behavior. The 
top pair of curves corresponds to the minimal g required for the massless case to achieve efficient 
preheating, p x ps p<p. As predicted by the heuristic picture, increasing m x or /(x) or decreasing g 
greatly suppresses amplification. 
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m x / m 

FIG. 3. A plot of g vs. ra x jm^ indicating the boundary between efficient preheating and 
inefficient preheating. Numerical results (circles) and the prediction based on the Mathieu equation 
(solid curve) are in good agreement. Note the lower bound on g required to have efficient preheating, 
a bound which increases as m x increases. 
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FIG. 4. For simplicity, the calculations in all three panels of the figure are done for a single 
k = mode of the x held. Figure (a) shows the number of %-particles (n x ) produced as a function 
of time for a typical case of barely efficient preheating. Each time 4> oscillates through its minimum, 
there is a sharp spike; however, particle production is measured by the time-averaged increase in 
n x . Initially, when q is large (t/2ir < 8), the increases in n x occur in discrete jumps from one side 
of a spike to another. This behavior is just like the static Mathieu equation for high q, as shown in 
(b). At later times (t/27r > 8), the growth is continuous between spikes and across them, similar 
to the static Mathieu equation at small q, as shown in (c). For t/2-K > 12, there are also rises and 
plateaus stretching many oscillations as q spends many oscillations in one instability regime. In 
all cases, the figures show that the growth has an exponential envelope with exponent v reliably 
predicted by the Mathieu picture. 
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